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Analysing panel flutter in supersonic 
flow by Hopf bifurcation 


Z. Monfared* and Z. Dadi 


Abstract 


This paper is devoted to the study of a partial differential equation (PDE) 
governing panel motion in supersonic flow. This PDE can be transformed to 
an ODE by means of a Galerkin method. Here by using a criterion which is 
closely related to the Routh-Hurwitz criterion, we investigate the mentioned 
transformed ODE from Hopf bifurcation point of view. In fact we obtain a 
region for existence of simple Hopf bifurcation for it. With the aid of Matlab 
and Hopf bifurcation tool, flutter and limit cycle oscillations of panel are 
verified. Moreover, Hopf bifurcation theory is used to analyse the flutter 
speed of the system. 


Keywords: Panel flutter; Limit cycle; Hopf bifurcation; Routh-Hurwitz cri- 
terion; Vibrations. 


1 Introduction 


Aerodynamics is a branch of dynamics concerned with studying the motion 
of air, particularly when it interacts with a solid object such as an Aircraft 
structure. 

On the other hand, flow-induced structural vibration is one of the most 
technical problems affecting the reliability, cost and safety of aircraft struc- 
tures. The vibration caused by a fluid flowing around a body is known as 
flow-induced vibration. Flow-induced vibrations best describe the interac- 
tion that occurs between the fluid’s dynamic forces and a structure’s inertial, 
damping, and elastic forces. The study of flow-induced vibrations has rapidly 
developed in aeronautical and nonaeronautical engineering. In aeronautics, 
flow-induced vibration is often referred to as flutter. Flutter is the instability 
of aeronautics structures under unsteady aerodynamic loadings. Panel flutter 
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is a phenomenon of self-exciting vibrations of skin panels of flight vehicle at 
high flight speeds resulting from the interaction between an elastic structure 
and the flow around the structure. Such vibrations typically have high ampli- 
tude and cause fatigue damage of skin panels. This flutter phenomenon was 
first observed during World War II; however, formal studies did not appear 
until the 1950s. Supersonic panel flutter is a key design consideration for 
some high-speed aerospace vehicles like spacecrafts and missiles, [14, 17, 20]. 
Moreover, for nonlinear systems flutter is usually interpreted as a limit cycle 
oscillation (LCO). If a limit cycle is created the system will oscillate forever. 


Aeroelastic flutter is a catastrophic structural failure, which needs to be 
avoided within the flight envelope of an aircraft structure. Hence, engineering 
researchers have paid much attention to studying the flutter and limit cycle 
motions of thin panels in recent years (see [2, 12, 14, 15, 16, 17, 19, 20] ). 

Furthermore, Hopf bifurcation theory can be utilized as an important tool 
for the determination of the flutter and limit cycle vibrations of panels. In 
addition, Hopf bifurcation tool can be used to analyze the flutter speed of 
the system. Hence, with the use of thin panels in shuttles and large space 
stations, nonlinear dynamics, bifurcations, and the chaos of thin panels have 
become more and more important. In the past decade, researchers have made 
a number of studies into nonlinear oscillations, bifurcations, and the chaos of 
thin panels. Holmes [5] studied flow-induced oscillations and bifurcations of 
thin panels and gave a finite-dimensional analysis. Then based on the analysis 
in [5], Holmes and Marsden [7] considered an infinite-dimensional analysis for 
flow-induced oscillations and pitchfork and fold bifurcations of thin panels. 
Holmes [6] then simplified this problem to a two-degrees-of-freedom nonlinear 
system and used center manifolds and the theory of normal forms to study 
the degenerate bifurcations. Yang and Sethna [18] used an averaging method 
to study the local and global bifurcations in parametrically excited, nearly 
square plates. From the von Karman equation, they simplified this system 
to a parametrically excited two-degrees-of-freedom nonlinear oscillators and 
analyzed the global behaviour of the averaged equations. Based on the stud- 
ies in [18], Feng and Sethna [3] made use of the global perturbation method 
developed by Kovacic and Wiggins [8] to study further the global bifurcations 
and chaotic dynamics of a thin panel under parametric excitation, and ob- 
tained the conditions in which Silnikov-type homoclinic orbits and chaos can 
occur. Zhang et al. [19] investigated both the local and global bifurcations 
of a simply supported at the fore-edge, rectangular thin plate subjected to 
transversal and in-plane excitations simultaneously. 

In this paper, a problem of flow-induced oscillations, that of panel flutter 
is considered. In fact here we investigate a partial differential equation which 
describes panel motion and obtain a region for the existence of a special type 
of Hopf bifurcation for it. This type of Hopf bifurcation occurs where a pair 
of complex conjugate eigenvalues of the Jacobian matrix passes through the 
imaginary axis while all other eigenvalues have negative real parts. Further- 
more, the existence of this type of Hopf bifurcation leads to flutter and limit 
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cycle motions of the panel which can cause failure of the structure. To the 
best of our knowledge, it is the first time that such a region for the existence 
of periodic solutions and Hopf bifurcation is being investigated. Moreover, 
by means of Matlab and the fourth and fifth-order Runge-Kutta (RK-45) 
method we do some numerical simulations. These simulations present our 
theoretical results and flutter and limit cycle oscillations of thin panel. More- 
over, the flutter speed is obtained by using Hopf bifurcation tool. 


2 Preliminaries 


In this section, we state some mathematical concepts and basic results. 


2.1 A criterion based on the Routh-Hurwitz criterion 
for the existence of simple Hopf bifurcation 


Important criterion that gives necessary and sufficient conditions for all of 
the roots of the characteristic polynomial to lie in the left half of the complex 
plane is known as the Routh-Hurwitz Criterion. This criterion is stated in 
the next theorem, see [16]. 


Theorem 2.1. (Routh-Hurwitz Criterion). Consider a polynomial of the form 


k— 


k 1 k-2 
Anz” + Ap_12 + Ap—22 + .. +40, 


the roots of this polynomial lie in the open left half-plane if and only if all the 
leading principal minors of the k x k matrix 


ay ao 0 
a3 a2 0 

Q= : ; 
G2k—1 @2k—-2 *** Gk 


are positive and ay > 0; we assume that aj =0 if 7 <0 orj >k. 


In dynamical systems, a bifurcation occurs when a small smooth change 
made to the parameter values (the bifurcation parameters) of a system causes 
a sudden qualitative or topological change in its behaviour. In general, at a 
bifurcation point, the local stability properties of equilibria, periodic orbits or 
other invariant sets change. Moreover, in many applications we are concerned 
with a special type of Hopf bifurcations, where a pair of complex conjugate 
eigenvalues of the Jacobian matrix passes through the imaginary axis while 
all other eigenvalues have negative real parts. These are called simple Hopf 


4 Z. Monfared and Z. Dadi 


bifurcations in [10], in order to distinguish them from the Hopf bifurcations 
with some other eigenvalues on the right half plane. Now consider the system 


X=f(X,y), X eR", per. (1) 


We use the notation in [4, 13] and mention the following theorem which states 
the sufficient conditions for existence of Hopf bifurcation. 


Theorem 2.2. Suppose that system (1), has an equilibrium (ao, fo) at which 
the following properties are satisfied: 


(SH1) Dz fpo(Lo) has a simple pair of pure imaginary eigenvalues and other 
eigenvalues have negative real parts. Therefore, there exists a smooth 
curve of equilibria (ax(1), 4) with (49) = xo. The eigenvalues 
Mw), ACM) Of Defuo(a(W)) which are imaginary at w = po vary 
smoothly with uw. Furthermore, if, 


(SH2) £(ReAw)| =a 40. 


H=Ho 


then Hopf bifurcation will occur. 


Proof. See [4]. 


Even though numerical computation of eigenvalues is feasible, it is ideal 
to have a criterion stated in terms of the coefficient of the characteristic poly- 
nomials rather than the traditional Hopf bifurcation criterion which is based 
on the property of eigenvalues. Specially for higher dimensional systems with 
many parameters, this criterion will be more convenient. See [10]. 

We denote the characteristic polynomial of the Jacobian matrix J() of (1) 
as: 
P(A; pb) = det(AIn — I(t) = pol) + pi(@aA +. + Pr(e)r”, 


where every p;(j) is a smooth function of u, and p,(j) = 1. And we consider 
the case po(ps) > 0, because there is not any nonnegative real root. Let 


pi(p) pol) + 0 
p3() po(w) -- 0 
Ln 7) = : . os . ’ 


P2n—1(H4) P2n—2(H) e Pn (KW) 


where p;(u) = 0 if ¢ < 0 or i > n. Moreover, when po(js) > 0 by the R-H 
criterion the polynomial P(); 1) of X has all roots with negative real parts if 
and only if the following n principal sub determinants of L,,(j1) are positive: 


e Di(u) = det(Li(u)) = pi(u) > 0 


pil) pole) 


© D2(y) = det(Lo(u)) = det ae pa() 


)>0 
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e Dr(u) = det(Ly(p)) > 0. Since Dp (i) = pn()Dn—i() and in our case 
Pn() = 1 the R-H criterion conditions can be stated as 


e Do( Lt) > 0, D,>0, Dz > 0, a, Dy_1 > 0. 


Theorem 2.3. Assume there is a smooth curve of equilibria (x(j), u) with 
x(Uo) = Lo for (1). Then conditions (SH1) and (SH2) for a simple Hopf 
bifurcation are equivalent to the following conditions on the coefficients of 
the characteristic polynomial P(A; 1): 


(1) po(uo) > 9, Di(fo) >0, «-, Drn—2(uo) > 9,  Dn—-1(so) = 0 


(2) (Pee te) #0. 


Proof. see [10]. 


2.2 Galérkin method 


Suppose we wish to solve the following boundary value problem of partial 
differential equations over the interval a < z < }, 


Lly(z, t)] + f(z,t) = 0, 


y(a,t) = 2a, 4(b,t) = Zp. 


A Galérkin method is used to approximate the problem by a sequence of 
finite dimensional problems. In other words, we consider the problem as a 
flow defined on a space V and then choose the finite dimensional subspace 
Vy C V of dimensional N and project our problem onto Vy. Reducing 
the problem to a finite dimensional vector subspace allows us to numerically 
compute uy (the solution of PDE) as a finite linear combination of the basis 
vectors in Vy. 

Now, let {y; an be an orthonormal basis of the finite dimensional sub- 
space Vy that satisfy the boundary conditions of the problem. Therefore, we 
can write 


N 
un(z,t) = Dd a;(t)p;. 


Through the use of orthogonal functions the error function Ey, represent- 
ing the difference between the exact and approximate solution, is minimized 
such that 


b 
/ Ey (z,t)p; dz=0, Vj =1,2,...,.N, 


where, 
En(z,t) = L[un(z,t)] + f(z, 6). 
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Each term in above equation gives an ODE in time for the N coefficients 
{a;(t)}fu4 and these must be solved numerically. For more information see 


(5, 7]. 


2.3 Dynamic pressure 


In fluid dynamics, dynamic pressure (indicated with g, or Q, and sometimes 
called velocity pressure) is the quantity of air measured by most airspeed 
instruments and defined by 


ere 
T= 5"; 


where p and v are density and velocity of the flow respectively. See [1], 
Section 3.5. 


3 Formulation of the problem 


Consider a supersonic stream of fluid passes above a thin plate with the 
length 1, fixed at the edges z = 0 and z = 1 .The panel is simultaneously 
subjected to an in-plane tensile load [. The fluid velocity is characterized 
in terms of the dynamic pressure qg, see Figure 1. Using nondimensional 
quantities, and assuming that the panel bends in a cylindrical mode (so 
that w(z,y,t) = w(z,t) is independent of y), the following nonlinear partial- 
differential equation, which is essentially a one dimensional version of the von 
Karman equations is considered for a thin plate: 


1 1 
Wee + AWtzz22 + gow; — {+ k fo wedz+ o fo We Wizdz} We, (2) 
+Weeez + QW, = 0, 


see [4, 7]. Here w = w(z,t) is the transverse displacement of the panel, 
a,o > 0 are (linear) viscoelastic damping parameters associated with the 
panel, 6 > 0 represents fluid damping and k& > 0 is a measure of the nonlinear 
axial (membrance) restoring forces generated in the panel due to transverse 
displacement. Moreover, all these parameters are assumed to be fixed, except 
q which can vary. 
Now, consider equation (2) with the following simply supported boundary 
conditions 

w(0,t) = wzz(0,t) = w(1,t) = wz (1,t) = 0. (3) 


A plate subjected to a compressive in-plane load with fluid flow over its sur- 
face may undergo complex motions resulting in dynamic instabilities (flutter) 
and associated limit cycle motions. 
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Figure 1: The panel flutter problem 


4 Bifurcation analysis 


It does not seem possible to solve the equation (2) explicitly. As it is men- 
tioned in [7] by using a Galérkin method the partial differential equation (2) 
together with the boundary conditions (3) can be transformed to a ordinary 
differential equation. In fact in [7] because of the simply supported boundary 
condition (3), the following family of orthogonal basis 


{9;(2)}j1 = {sin(jnz)} Ly 


is chosen. Then by writing , w(z,t) = > 3 a,(t)y;(z) and applying the 
Galérkin procedure for N = 2 (two modes) to the governing Equation (2) 
and using the orthonormality of the bases and the relationships 


1 1 
Wi —_— hi MW 
| Wy -Ws a= [ ww, dz, 
0 0 
1 1 
" _ roy 
| W; Ws dz = -{ wW;-W, dz, 
0 0 


the following ODE in the time dependent amplitude coefficients a;(t) is ob- 
tained. 


=A rt+f(z), zcéR*t,qeR, (4) 
where 
ay 
a2 
c=| “|, 
ay 


a2 
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0 0 1 0 
Br 0 0 0 1 
q | w?7(T — 1) 8q/3 —(an* + \/q6) 0 : 
—8q/3 4n?([ — 4x?) 0 —(16a7* + \/96) 

and 

0 

0 

f a, —fi ’ 
2 

in which 


4 : ; 
iy= 5 {k(aq 4a3) + o(a1a1 + 4a2d2)}az, 
fr = 2n*{k(az + 4a3) + a(a1a1 + dagd2)}ao. 


Moreover, a; = x;(i = 1,2) are the amplitudes of normal two modes. 

Now we investigate Hopf bifurcations of system (2) for the trivial equilib- 
rium position z = 0 or w(z,t) = 0 by using a criterion which is closely related 
to the Routh-Hurwitz (R-H) criterion to obtain a region of simple Hopf bifur- 
cation parameters. For this purpose, we mention the sufficient conditions for 
existence of simple Hopf bifurcation of system (4) in the following theorem. 


Theorem 4.1. Suppose that for q = qo and the trivial equilibrium position 
x =0 of (4), the following relations satisfy: 


1) po(qo) > 0 


2) pi(qo) > 0 


P1(qo) Po(qo) 
3) det >0 


p3(qo) P2(qo) 
Pi(qo) Po(go) 0 


4) D3(qo) = det | p3(qo) p2(qo) Pi(qo) | = 0 


0 1 ps(qo) 
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po(qo) = 4n2(U ~ n2)(C — 4?) + S20, 

pi(qo) = oN n*)(16a74 + /q6) + 4(P — Ar?) (an* + /09)}, 
p2(do) = (an* + \/G95)(16an* + 7d) — 12(5T — 1727), 

p3(qo) = 17ar* + 2,/q06. 


Then qo is a simple Hopf bifurcation value for system (4) at the trivial equi- 
librium position x = 0. 


Proof. By computing the characteristic polynomial of (4) and using theorem 
2.3., the above assertion can be proved. 


The panel transverse displacement and velocity for two modes (N = 2) 
can be evaluated by the following equations 


t)= ay (t)sin(j7z), (5) 
= x a, (t)sin(jrz). (6) 


Therefore, by the above theorem we can find a region for existence of simple 
Hopf bifurcation for equation (4) and therefore for equation (2). 


5 Numerical simulation 


Here numerical simulations are carried out to support our theoretical results 
and show panel flutter. 


Example. Consider system (4), by the aid of software Auto for 6 = 1,T = 
2,a = 0 and q ~& 323.24789021 we can show that this system undergoes 
Hopf bifurcation at the equilibrium x = (21, 2%2,%73,x%4) = (0,0,0,0). More- 
over, these values of parameters by [4, 5, 7] are physically meaningful. In 
addition, for the mentioned values of parameters the characteristic polyno- 
mial of the Jacobian matrix of (4) at x for qo is: 


P(A; qo) © 857953.77463 + 27998.10001A + 1880.50639A2 + 35.95819\3 + \4. 


So 
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Figure 2: The related trajectories in the phase space 21 —x2 and time history for different 


values of q. Initial condition is very close to the equilibrium x = 0: (a)q = 300; (b) ¢ = 327; 


(c) q = 1000 
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Figure 3: Phase Portrait in the space w — wt and time History of the panel transverse 
displacement for different values of g : (a)q = 300; (b) g = 327; (c) g = 1000 
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po (qo) © 857953.774634701 > 0, 
D1 (qo) = p1(qo) & 27998.1000129 > 0, 


p1(qo) po(qo) 


Data) =a t( 
2(qo) . p3(qo) p2(qo) 


) ~ 21800139.664474831 > 0, 


9 


Pi(qo) Po(qo) 0 
Ds3(qo) = det (mio p2(qo) nin 
0 1 ps(qo) 


and 


dD 
a ~ —4937513.3537057098 # 0. 


Therefore, by Theorem 4.1., gg ~ 323.24789021876 is a simple Hopf bifur- 
cation value for (4) at 2. Furthermore, for go the Jacobian matrix Ag at x 
has a pair of pure imaginary eigenvalues Az, A2 ~ £27.90392896678841892, 
and a couple of complex conjugate eigenvalues with negative real part 
Rer3 = RerA4, ~ —17.979095923287. 

We solved the equation (4) by means of the fourth and fifth-order Runge- 
Kutta (RK-45) method. Then by numerical simulations we showed the occur- 
rence of Hopf bifurcation. The related trajectories in the phase space x1 — £2 
(the space of the amplitudes of two modes) and the time history responses 
for different values of parameter q are presented in Figure 2. 

Moreover, by equations (5) and (6), the panel transverse displacement 
(w(z,t)) and velocity (w,(z,t)) for two modes can be determined. So, the 
related trajectories in the phase space w — w; and the time history responses 
of the panel transverse displacement for different values of parameter q are 
illustrated in Figure 3. 

As it is illustrated in Figures 2 and 3, the equilibrium point x = 0 or 
w(z,t) = 0 is a stable focus when the dynamic pressure of flow is less than 
qo. In this case by passing the time the amplitude of panel vibrations will 
vanish. While for g > qo the equilibrium point x = 0 or w(z,t) = 0 turns 
out to be an unstable focus surrounded by a stable limit cycle. That means 
that by passing the time, the amplitude of panel oscillations will increase and 
finally panel will vibrate with a fixed period for ever. Furthermore, in this 
case for large values of the dynamic pressure q, panel vibrations can have 
high amplitude and cause catastrophic failure of the structure. This increase 
of amplitude will be faster for larger values of q. 

In addition, since the Hopf bifurcation occurs for the dynamic pressure 
do © 323.24789021876 Pa, due to Section 2.3., by knowing the density of flow 
(air) the flutter speed can be obtained. For example at sea level and at 15 
°C’, air has a density of approximately 1.225 Fg. Hence, for this density of 
air the flutter speed is approximately 22.97284609 in the absence of linear 
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viscoelastic damping parameter a. In Figure 3 our simulations show flutter 
and limit cycle oscillations (LCO) of the panel without considering the vis- 
coelastic damping parameter a. 


6 Conclusion 


In this paper, we extended the previous results in [5, 6, 7, 8, 17, 18, 19] to 
study the vibrations of a thin panel fixed at two edges. Indeed we found a re- 
gion for existence of simple Hopf bifurcation for a partial differential equation 
governing panel motion. Because, the existence of simple Hopf bifurcation 
can lead to flutter and limit cycle oscillations of the panel. Numerical sim- 
ulations were carried out by using the fourth and fifth-order Runge-Kutta 
method, to support our analytical results. In fact by simulations and Hopf 
bifurcation theory, we showed the occurrence of flutter and limit cycle mo- 
tions of thin panel. Then Hopf bifurcation tool was used to calculate the 
flutter speed of the system. Moreover, numerical simulations presented vi- 
brations of thin panel can have high amplitude which cause damage of panel. 
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